Dynamical response functions and collective modes of bilayer graphene 
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Bilayer graphene (BLG) has recently attracted a great deal of attention because of its electrically 
tunable energy gaps and its unusual electronic structure. In this Letter we present analytical and 
semi-analytical expressions, based on the four-band continuum model, for the layer-sum and layer- 
difference density response functions of neutral and doped BLG. These results demonstrate that 
BLG density-fluctuations can exhibit either single-component massive- chiral character or standard 
two-layer character, depending on energy and doping. 
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Introduction — Recent progress [1 in the isolation and 
experimental exploration of large area single and mul- 
tilayer graphene systems has opened up a new topic in 
two-dimensional electron system (2DES) physics. These 
atomically thin 2DESs exhibit a rich variety of unique 
properties that are presently under active investiga- 
tion. In particular, the peculiarities of one (SLG), two 
(BLG) [21 131 SI |5| , and three layer systems are quite dis- 
tinct. This Letter is motivated by on-going experimental 
work on suspended BLG [6 and improved BLG samples 
on SiC [7]- We anticipate that electron-electron interac- 
tions will have a crucial influence on the emerging Fermi 
liquid, collective excitation, angle-resolved photoemission 
spectroscopy (ARPES) [7], and tunneling properties of 
BLG systems. 

Many-body effects in BLG have been studied by sev- 
eral authors [H [H [lOl [Til H H [H] • However density- 
response functions, which are the starting point for 
detailed many-body physics considerations in charged- 
particle systems, have so far been calculated [12] only in 
the static limit, only in the density-density channel (see 
below), and only within a two-band model [15] whose 
applicability is limited to low-densities and low-energies. 
Dynamical screening and collective effects in BLG are 
still largely unexplored. 

In this Letter we present analytical expressions based 
on the full four-band continuum model for the dynami- 
cal susceptibilities of undoped BLG and semi-analytical 
expressions for the same quantities in doped BLG. Our 
results provide the key technical ingredient necessary 
for many-body theory calculations that are based on 
the random-phase-approximation (RPA) or its general- 
izations, whether directed toward thermodynamic quan- 
tities (like charge and spin susceptibilities) or toward 
quasiparticle dynamics [21 17] . They exhibit many inter- 
esting features which foreshadow key aspects of many- 
body correlation physics in these systems. In the present 
paper we present detailed RPA predictions for the col- 
lective plasmon excitations of BLG, which are expected 
to be directly observable in electron energy loss spec- 



troscopy studies and, as in the SLG case, are responsible 
for the many-body features observable in ARPES spec- 
tra [HHi]. 

Four-hand linear-response theory — BLG is modeled as 
two SLG systems separated by a distance d and cou- 
pled by both inter-layer hopping and Coulomb inter- 
actions. Most of the properties we discuss below de- 
pend qualitatively on the Bernal stacking arrangement 
in which one sublattice (say A) of the top layer is a near- 
neighbor of the opposite sublattice (say B) of the bottom 
layer. Neglecting trigonal warping, which is important 
only at extremely low densities and is presently masked 
by uncontrolled disorder, the single-particle Hamiltonian 
is (n = 1) r = Efc,a,/3 4,a^«/3(^)cfc,/3, where T{k) = 
-V7^7^7 • k - t±{j^j^ + i^y)/2. Here v (~ 10^ m/s) 
is the Fermi velocity of an isolated graphene layer, 
(~ 0.3 eV) is the inter-layer hopping amplitude, and the 
7^ are 4x4 Dirac 7 matrices in the chiral representa- 
tion [18^ (7^ = —^7^7"^ 7^7^). The Greek indices ac- 
count for the sublattice degrees of freedom in top (1 = A, 
2 = B) and bottom (3 = A, 4 = 5) layers. Two elec- 
trons in the same (S) layer interact via the 2D Coulomb 
potential Vs(^) = 27re^/(eg). Electrons in different (D) 
layers interact via Fd(^) = ^s(^) exp {—qd). 

For response function calculations it is convenient to 
work in the single-particle Hamiltonian eigenstate basis. 
Diagonalization of T{k) yields four hyperbolic bands [8] 
(see Fig.[l]) with dispersions, £1,2 (^) = + ^1/4+ 

t±/2 and €s,4{k) = ±^/v^k^ ^t\/4 - t±/2. In this ba- 
sis the interaction contribution to the Hamiltonian is 
Hint = {^S)-'E^[V^{q)p^p_^ ^ V_{q)f^f_^], where 
S is the 2DES area, V± = {Vs ± Fd)/2, and pq 
and Tq are respectively the operators for the sum 
and difference of the individual layer densities: Pq = 

T.k,X,y ^k-q,xi^k-q,k)xX'Ck,X' with = 

and fq = 5])fc,A,A' cl_q y^{Sk-q^k)xX'Ck,X' with Sk-q,k = 

Ul^_^^^Uk' Here Uk is the unitary transformation from 

sublattice to band labels A, A' [8]. From TYint we thus see 
that two response functions are necessary for the eval- 
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FIG. 1: (Color online) BLG continuum- mo del band structure. 
Depending on the doping level, the 2DES can have either one 
or two conduction band Fermi surfaces. 



uation of collective modes and ground-state properties 
of BLG: the total-density response function, Xpp{q^^) = 
{{pq; p-q))uj/S, and the density-difference response func- 
tion Xtt(<7,^) = {{fq;f_q))^/S. Here {{A;B))^ is the 
Kubo product [19, 20 . 

Noninteracting response functions and RPA screening — 



In the noninteracting limit the linear-response func- 
tions introduced above have the standard eigenstate- 
representation form [19j: 



A, A' ' 



(27r)2 Z + Afc,A;fc',A' 



Mh,\-k'^' , (1) 



where z = uo ^ zO+, k' = k -\- n^^x are band occu- 
pation factors, and Aj^^x;k',x' = £k,x — ^k' ,x' are band- 
energy differences. Here A^fc,A;fc',A' is \{Vk^k>)xx'\^ for the 
total-density response and \{Sk^k')xx'\^ for the density- 
difference response. We have evaluated X^p^TT) ~^ 
/^pp(TT) analytically for undoped BLG, i.e. for the case 
in which bands 2 and 4 are full and bands 1 and 3 are 
empty. Here we report only results for the imaginary 
parts of these response functions. The corresponding 
analytical expressions for the real parts, which can be 
derived from a standard Kramers-Kronig analysis, are 
extremely cumbersome and will be presented elsewhere. 
After very lengthy algebra we have reached the following 
results (per spin and per valley): 
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y , +2V^((7,cj,Cc;) \g{q,u;,uj-)\ 



e{g{q,u;,u;) - tl_) 



uj^g{q,uj_,uj_) - \g{q,uj,uj_)\ Q{g{q,uj_,iJ_) - t\) \ ^ I ... 



(2) 



and 
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v^P{q,u_)-2v^q^ -2t\ 



+ 2\/^(g,^-,cc;_) \g{q, uj, uj- 



uj^g{q,ij,ij) - \g{q,ij,uj_)\ Q{g{q,ij,ij) - t\)\ ^ { ... 



Q{g{q,uj-,uj-) - t\) 



(3) 



where a;± = ± t^, g{q-, ^) = — v'^q^ ^ f{q^ co) = 
q^y[g{q,^J,^J) —tj]/g{q^uo^uo)^ and 0(x) is the usual step 
function. Eqs. ([2| and ([3| greatly simplify the analysis of 
many-body effects in BLG and are an important result 
of this work. 



For both density-sum and density-difference channels, 
the response functions of the doped system can be written 
as x^^^ = x^^^^ + ^X^^^- We find that the corrections due 
to doping can be reduced to simple but cumbersome ID 
integrals: 



t±/i2v) 



. I [dppiTT) (5, z, -t±) + jpp(TT) {s,q,z,t±)]ds \ ^ < ... \ , (4) 



3 



where 



a{z)a{z+)]\gn{lRe[P{q,z)]) 



R{q) 



4[aHz+) 



(5) 



jpp = sgn{^e[P{q,Z-)]) 



^/Q{q,z_) R{q) 
a?{z-)-z'^ 4[a2(z_)- 



4(s - 1±/2) 



(6) 



and 



^TT = sgn(3?e[P((7,z)])- 



R{q) 



Jrr 



a?{z+)-z'^ A[a?{z+) - z"^] 4 4(s + t_L/2) 
-[v'^q^ + t^z-a?{z_)fsgn{^e[P{q,z_)]) R{q) 5 z 



\z.)-z^yQ{q,z.) 



4[a^{z-)-z^] 4 4{s-t_L/2) ' 



(7) 



(8) 
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FIG. 2: (Color online) BLG static response in units of 
the Fermi- level density-of-states of band 3, u = (sf + 
t±/2)/(27rv^), as a function of ^//cps- a) Xpp^(^, 0). The 
(black) dashed line is the result obtained within the two- 
band model [12 , while the (red) solid line is the result ob- 
tained within the four-band model for doping level n = 
10^^ cm~^, corresponding to the (red) solid line in Fig. [l] 
The (blue) dash-dotted line gives the static response for 
n = 5 X 10^^ cm~^ corresponding to the (blue) dash-dotted 
line in Fig. [l] Inset: small momenta region of the heavily- 
doped result. From left to right, the vertical dashed lines are 
at 2/cfi, k^i + A;f3, and 2/cf3. b) Xtt(^^^) with the same la- 
beling as in panel a). The two-band- model Xtt(^' 0) reported 
here has been calculated with a cut-off kc = t±/v. 



Here z± = z±tj_, a{z) = z-\-2vs, P{q^ z) = v'^q^ — za{z)^ 
Q{q,z) = v^q^ + a" {z) z^ + q^ [t\ -a\z)-z% and R{q) = 



\4v^q^ 



Note that the first term inside the 



curly brackets in Eq. Q is finite only if the high-energy 
band ei{k) is occupied {i.e. only if the Fermi energy 
Sf > t±). Eqs. (|4|-([8| constitute the second important 
result of this work. 

The static limit of these response functions 
/^pp(TT)(^'^ = 0) is illustrated in Fig. 2 for both 
lightly and heavily doped bilayers [21 . In the low- 



density limit Xpp{q^^) exhibits a strong Kohn anomaly 
at g = 2/cf3 associated with [12 the massive-chiral 
behavior of band 3 at energies below ~ tj_. We see 
in Fig. [2] that response functions calculated in the 
two-band model [12] (dashed line in Fig. [2]) overestimate 
the strength of this non-analyticity because they do not 
capture the gradual change in the single-particle eigen- 
state character of band 3 from the coherent two-layer 
wavefunctions at low energies to weakly coupled SLG 
wavefunctions at high energies. For the same reason the 
two-band model completely misrepresents the large-g' 
behavior, failing to capture the linear increase in x^^^ 
at large q which closely mimics SLG behavior. In the 
high-density limit Xp^p{Qi^) becomes rather similar to 
its SLG counterpart. The Kohn anomaly at 2/cfi, which 
still has BLG character at this energy, is relatively 
strong while the anomaly at 2/cf3, which already has 
more single-layer character, is strongly suppressed. 
The real-space Friedel oscillations (FOs) exhibit cor- 
responding changes [22] as the occupation of band 1 
increases at high densities. In panel b) we clearly see 
that the two-band model is even more inadequate in the 
density-difference channel. (In fact the integrals which 
appear in Xtt(^'^) ^^^^ ultraviolet divergence |^ 
in the two-band model.) Xtt larger than Xp^p at small 
q for low-densities because of the two-layer character 
wavefunctions are easily polarized. At higher densities 
Xyt Xp^p are nearly identical, as expected when the 
two-layers respond nearly independently. 
RPA theory of collective modes — The RPA response func- 
tions of the interacting doped system are given by 



Xpp(TT) 



v±x 



(0) 



(9) 



The interacting-system susceptibilities are determined by 
the density n, d (which we have taken tohe d = 3.35 A), 
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FIG. 3: (Color) BLG RPA dynamical dielectric functions for 
low-doping n = 10^^ cm~^ and aee = 0.5. The left panel 
shows [l/epp{q, cu)] as a function of q (in units of /cfs) and 
(jj (in units of £f). The right panel show [1/£tt(^, <^)]. 
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FIG. 4: (Color) High density (n = 5 x 10^^ cm"^) BLG RPA 
dynamical dielectric functions for the same model parameters 
as in Fig. [3] The inset in the right panel highlights the inter- 
subband plasmon. 



t± (which we have taken to be 0.35 eV), and by the ef- 
fective fine-structure constant o^ee = /{hve). In Figs, [s] 
and [4] we plot the imaginary parts, [1/ Spp^q^uj)] and 
[1/£tt(^, ^)], of the inverse dynamic dielectric func- 
tions which provide a portrait of BLG density-sum and 
density-difference fluctuations and collective modes. The 
collective fluctuation physics in the high-density limit is 
much like that of an ordinary bilayer [24 , as expected. 
A clear in-phase bilayer plasmon, whose frequency goes 
to zero like for g ^ appears at low-energies and 
is Landau-damped at relatively low frequencies by inter- 
band transitions. An out-of-phase inter-subband plas- 
mon appears in £tt(^,^) just above the transition fre- 
quency between bands 3 and 1. At low-densities, how- 
ever, these results show that the collective fluctuation 
physics in BLG is quite unusual. The inter-subband 
plasmon is Landau damped at all wavevectors for < 
t±/2^ i.e. for densities below ric = 3{t±/v)'^ /{Att) ~ 
7 X 10^^ cm~^, which is smaller than the critical den- 
sity ni = 2{t±/v)^/7r - 18 X 10^^ cm-^ at which the 
band €i{k) is populated. The density-sum plasmon still 
appears and still has dispersion but the physics which 
determines the coefficient of ^yq is completely altered [25 . 
The shark-fin structure around uo = and q — \f2 k^s 
inside the e-h continuum is a direct consequence of the 
J = 2 massive chiral fermion behavior because it leads to 
suppressed scattering from a state with momentum k to 
a state with final momentum {k-\-q) ± k. The disappear- 
ance of the inter-subband plasmon in [l/err{Q^^)] 
occurs because the mode would be Landau damped even 
if strong, and because long wavelength transition am- 
plitudes between bands 1 and 3 are suppressed at low- 
energies. 

In summary we have demonstrated that the density- 
sum and density-difference fluctuations in BLG crossover 
from those of an unusual massive- chiral single-layer sys- 
tem to those of a weakly coupled bilayer as carrier- 
density, wave vector, and energy increase. The analytic 



and semi-analytic results for RPA response functions ob- 
tained here will simplify efforts to understand the many- 
body physics of this unique 2DES. 
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